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1 Introduction 


The dynamical state of the ocean and atmosphere is taken to be a large- 
dimensional random vector in a range of large-scale computational applications, 
including data assimilation, ensemble prediction, sensitivity analysis, and pre- 
dictability studies. In each of these applications, numerical evolution of the 
covariance matrix of the random state plays a central role, because this matrix 
is used to quantify uncertainty in the state of the dynamical system. Since atmo- 
spheric and ocean dynamics are nonlinear, there is no closed evolution equation 
for the covariance matrix, nor for the mean state. Therefore approximate evo- 
lution equations must be used. 

This article studies theoretical properties of the evolution equations for the 
mean state and covariance matrix that arise in the second-moment closure ap- 
proximation (third- and higher-order moment discard). This approximation 
was introduced by EPSTEIN [1969] in an early effort to introduce a stochas- 
tic element into deterministic weather forecasting, and was studied further by 
FLEMING [1971a, b], EPSTEIN and PITCHER [1972], and PITCHER [1977], 
also in the context of atmospheric predictability. It has since fallen into disuse, 
with a simpler one being used in current large-scale applications. The theo- 
retical results of this article make a case that this approximation should be 
reconsidered for use in large-scale applications, however, because the second- 
moment closure equations possess a property of energetic consistency that the 
approximate equations now in common use do not possess. A number of prop- 
erties of solutions of the second-moment closure equations that result from this 
energetic consistency will be established. 

Suppose the dynamics of the state s £ l w are given by a system of nonlinear 
ordinary differential equations, 

ds 

— + f(s, t) = 0, (1) 

where t is time, f : S x T -» R w , 5 C R n is a state space appropriate for Eq. (1), 
and T = [to,T] is a closed time interval. The initial condition s to is taken to 
be a random state, s to £ S with probability one, so that the problem to be 
solved is the stochastic initial- value problem for Eq. (1). Technical assumptions 
on f are stated in Secs. 2 and 5. In addition, it will be assumed that the dy- 
namics are conservative, and a nonlinear transformation is introduced in Sec. 3 
to ensure that the total energy conserved by solutions of Eq. (1), stochastic or 
deterministic, is just E = A || s 1 1 2 , where || • || denotes the Euclidean norm on l w . 

A simple sufficient condition, which is natural for conservative dynamics, 
under which the stochastic initial-value problem for Eq. (1) is well-posed, is 
stated in Sec. 4. Under this condition, the solution of the stochastic initial- value 
problem defines a second-order stochastic process, that is, one that has a mean 
and covariance matrix at each time t in a closed time interval. Furthermore, it 
follows immediately from conservation of total energy E for this process that 



where s = s t £ S is the mean state of this process, P = P t £ M. NxN is the 
covariance matrix of the process, and tr P is the trace, or sum of the diagonal 
elements, of P. This means that the uncertainty in the random state s, as 
measured by the total variance V = trP, can increase (decrease) only as a 
result of extracting energy from (inserting energy into) the mean state s, with 
the change in total variance balanced exactly by twice the change in total energy 
^ ]|s|| 2 of the mean state. The mean state and covariance matrix cannot be 
calculated from Eq. (1) without approximation, however, unless f is linear in 
s, and one would like to develop approximate evolution equations for s and P 
that in the nonlinear case at least preserve this basic conservation property. 
A closed system of ordinary differential equations for the mean and covariance 
matrix whose solutions satisfy Eq. (2) will be said to be energetically consistent , 
after FLEMING [1971a, p. 872], 

The second-moment closure equations for approximate evolution of s and P 
are the closed, nonlinearly coupled differential equations 


ds 

dt 


+ i(s,t) + 

3 k 


<9 2 f(s, t) 
dsjdsk 


Pjk = 0 , 


( 3 ) 


- + F(s,f)P + PF T (s,f) = 0, (4) 

where Pjk is the (j, k) th element of the covariance matrix P whose evolution is 
described by Eq. (4), F = di /ds is the Jacobian matrix of f , and the superscript 
T denotes transposition. The evolution of P in Eq. (4) depends on that of the 
mean state s given by Eq. (3), through the dependence of the Jacobian matrix 
on the mean state, and the evolution of the mean state also depends on that 
of the covariance matrix, through the double-summation term in Eq. (3). If 
f(s ,f) is linear in s, then the second partial derivatives of f with respect to the 
state variables all vanish, so that the double-summation term vanishes. Hence 
this term is called the nonlinear coupling term. Equations (3) and (4) are to be 
solved together for initial conditions s to and P to , with s to being the mean of the 
random initial condition s to for Eq. (1) and P to being the covariance matrix of 

Sto ' 

Conditions under which this initial-value problem is well-posed on a closed 
time interval are given in Sec. 5, where a stochastic process having the solution 
(s, P) of the initial- value problem as its mean and covariance matrix is also 
defined. It is shown in Sec. 6 that the solution satisfies Eq. (2). Thus the 
second- moment closure equations are energetically consistent. For quadratically 
nonlinear f, energetic consistency of the second-moment closure equations was 
established by EPSTEIN [1969] and studied in detail by FLEMING [1971a], who 
also established energetic consistency of the third-moment closure equations for 
quadratically nonlinear f. The energetic consistency result established in the 
present article holds for general f. 

The derivation of Eq. (2) given in Sec. 6 for the second-moment closure 
equations shows that the exchange of energy between the mean state and the 
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stochastic perturbations, which balances exactly, occurs solely through the non- 
linear coupling term in Eq. (3) and the symmetric part F 5 = i(F + F T ) of the 
Jacobian matrix in Eq. (4). In other words, rewriting Eq. (4) as 

dP 

— + F a (s, t)P - PF a (s, t) + F s (s, t) P + PF s (s, t) = 0, (5) 

dt 

where F a = i(F — F T ) is the anti-symmetric (skew-symmetric) part of F, one 
has immediately that 

+ 2 tr F s (s, t)P = 0, 

and it is shown in Sec. 6 that Eq. (3) gives 


4 ff 

dt 


2 tr F s (s, t)P = 0, 


with contribution only from the nonlinear coupling term, not from the term 
f (s, t) in Eq. (3). Equation (2) then follows. 

In case f(s ,t) is linear in s, then not only does the nonlinear coupling term 
vanish, but for conservative dynamics F s vanishes as well, so that Eqs. (3) and 
(4) become simply 


ds 

— + f(s,£) = 0, 

(6) 

dP 

— + F°P - PF a = 0, 

dt 

(7) 


with F a independent of s by linearity. Thus the effect of nonlinearity, to second- 
moment closure, is to introduce the nonlinear coupling term and the terms 
F s (s, t)P + PF s (s, 1), and these terms together are in energetic balance. Non- 
linearity also introduces dependence of F“ on the mean state, but Eq. (7) is 
energetically neutral, satisfying dtrP/dt = 0, regardless of any dependence of 
F a on the mean state. 

The approximation now widely used in large-scale atmospheric and oceanic 
applications is to retain Eq. (4) as it stands, but to neglect the nonlinear coupling 
term. This is the approximation made for instance in four-dimensional varia- 
tional data assimilation (TALAGRAND and COURTIER [1987]; COURTIER 
and TALAGRAND [1987]; THEPAUT et al. [1996]) and in a variety of algo- 
rithms based on singular vector calculations (e.g. BUIZZA and PALMER [1995]; 
MOLTENI et al. [1996]; MOORE and KLEEMAN [1997]). This approximation 
is convenient for computations, because the mean state can then be evolved in- 
dependently of the covariance matrix. Neglecting the nonlinear coupling term, 
however, destroys energetic consistency. Furthermore, the nonlinear coupling 
term and the terms in the covariance evolution equation all have formally the 
same order of magnitude, since all are linear in the covariance matrix. More- 
over, in the sense of contribution to the total variance, the terms retained in the 
covariance evolution equation that arise from nonlinearity, F S P+PF S , have pre- 
cisely the same magnitude, and opposite sign, as that of the nonlinear coupling 
term which is neglected. 
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The role of F s in the energetic coupling of the second-moment closure equa- 
tions is studied further in Secs. 7 and 8. It is shown in Sec. 8 that if f is genuinely 
nonlinear, as defined there, then F s has at least one positive and one negative 
eigenvalue. This means that when the dynamics are genuinely nonlinear, there 
is always a direction in state space in which uncertainty decays, as well as a 
direction in which uncertainty grows. One implication is that if f is genuinely 
nonlinear, then neglect of the nonlinear coupling term can lead either to increase 
or decrease of the perceived uncertainty. 

When the nonlinear coupling term is neglected, Eq. (4) is a linear equation 
to be solved once the mean state has been calculated, and so its solution can be 
expressed in the form 

P t = M t , to P to Ml to , (8) 

where M f to is the fundamental matrix (alternatively, solution operator, or tan- 
gent linear propagator) of the perturbation dynamics corresponding to Eq. (1) 
linearized about the mean state. This expression is particularly convenient for 
singular-value calculations in large-scale applications. 

When the nonlinear coupling term is retained, the solution of Eq. (4) can 
still be expressed in the form (8), but for a matrix M t>to that itself depends 
on the covariance matrix. In Sec. 9, Eq. (2) is used to establish simple time- 
independent upper bounds for trP t /trP 4o , which hold also for the covariance 
matrix of the original stochastic process. When the nonlinear coupling term 
is neglected, the largest singular value of M^ to is the least upper bound for 
trP t /trPt 0 , but in general there is no time- independent upper bound since 
Eq. (2) does not hold. 

A minimum requirement for solutions of Eqs. (3) and (4) to approximate the 
mean and covariance matrix of solutions of the stochastic initial- value problem 
for Eq. (1) is for s to and P to to be the mean and covariance matrix of some 
random state, s io £ S with probability one. Because ocean and atmospheric 
dynamics contain state variables that are constrained to be positive, such as 
mass and temperature variables, or to satisfy other constraints, this minimum 
requirement implies that P to cannot be chosen independently of s to in general. 
Thus the initial-value problem for Eqs. (3) and (4) must be posed with some 
care. Toward this end, the deterministic initial-value problem for Eq. (1) is 
reviewed briefly in Sec. 2, and the stochastic initial-value problem for Eq. (1) is 
described in some detail in Sec. 4. The initial-value problem for Eqs. (3) and (4) 
is posed in Sec. 5, and restrictions on the initial covariance matrix are illustrated 
there using a spatially discretized version of the one-dimensional shallow-water 
equations. 

Brief concluding remarks are given in Sec. 10. 
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2 The Deterministic Initial- Value Problem 


Suppose that the evolution of a real vector s £ 1 N is governed by a nonlinear 
system of ordinary differential equations, 

^+f(s) = 0, (9) 

where f(s) = f(s,t) may depend explicitly on time t. Suppose also that f : 
1 5 x T -» Ht N , where S C is a convex open set, that is, an open set in 
R^ such that the line segment between any two points in S lies entirely in 
S, and where T = [fo,Z] is a fixed, closed time interval; T — to is finite but 
may be arbitrarily large. The open set S is called the state space, an element 
sGSis called a state, or state vector, and the N components of a state are the 
state variables. If the actual system under consideration is complex, Eq. (9) is 
obtained by separation into real and imaginary parts. 

This section provides a brief review of the deterministic initial- value problem 
for Eq. (9), in which one is supposed to find, for each initial state s to £ S, a 
solution s = s(t) e S that satisfies s(f 0 ) = s to . Sufficiently strong hypotheses 
will be imposed on f to guarantee existence and uniqueness of solutions on a 
(possibly short) half-open time interval Z* = [fo,Z*) C Z, with Z* depending 
in general on s to , and also to guarantee that the solution is in the space C 1 ( %) 
of functions with one continuous time derivative on 7 *. The corresponding 
stochastic initial-value problem, in which s to will depend on a probability vari- 
able, is considered in Sec. 4. It should be noted here that Assumption 2 below is 
stronger than necessary for establishing uniqueness of solutions of Eq. (9), but 
that this assumption will be required later, in Sec. 5, to ensure that solutions 
of the covariance evolution equation have one continuous time derivative. 


Assumption 1 f £ C(S x T), the space of continuous functions on S x T. 


Assumption 2 F = F(s,t) G C(S x T), where ¥ = di/ds denotes the N x N 
Jacobian matrix of{, whose ( j,k) th element is given by 


F jk (s,t.) 


d/j(M) 

ds k 


That is, Fjk £ C(S xT) for j, k = 1, . . . , N. 


Assumption 1 guarantees that for each s to £ S, there exists a time interval 
% = [fo,Z*), with T* = T*( s to ) < T, and at least one solution s (t) £ S for all 
t £%, such that s(f 0 ) = s to . It also guarantees that every such solution is in 
C ,1 ( %). Furthermore, existence of a solution ceases only if it hits the boundary 
of 5 (e.g. CODDINGTON and LEVINSON [1955, Ch. 1, Thm. 4.1]), where f 
may not even be defined. Assumption 2, along with the fact that S was taken to 
be convex, implies that f satisfies on S x T a Lipsclritz condition in s, uniformly 
in t , and this in turn guarantees that there is at most one solution on any time 
interval. 
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In particular, if the state space S is all of K w , then S has no boundary, and 
so for each s io G S there exists a unique solution s(f) G S over the full time 
interval T, with s(fo) = s to , and this solution is in C 1 (T). More generally, the 
same is true if the state space is not all of 1^, that is, if Assumptions 1 and 2 
hold only for a proper subset S of R N , provided the dynamics (9) are such that, 
for some choice of the state space, no solution can ever hit the boundary of S 
by time T. It is often the case that the physical problem at hand dictates that 
the state space cannot be all of R w , but that there is a choice of state space for 
which all solutions exist in S over the full time interval T, as discussed further 
in Sec. 3. 

Hereafter, the unique solution s(f) G S of Eq. (9) on some time interval 
% = %(s to ) C T, such that s (to) = s to , or on all of T in case it exists for all 
time t, G T, will be denoted by s t . The continuous path through state space 
traced by s t as time progresses is called the trajectory corresponding to s to . The 
bold letter s without a subscript will usually denote an arbitrary point in the 
state space, or in 1^. 

An implication of Assumption 2 beyond uniqueness of solutions, which is 
important for application to the covariance evolution equation, is that solutions 
depend continuously on parameters such as initial conditions. Regarding each 
trajectory s t as a function of its initial point s to , one finds that the TV x N 
matrix M = M(f) = ds t /ds to , whose (j, k) th element is given by 


Mjk = 


d(s to )fe’ 


satisfies the simple linear equation 


— b F(st, £)M = 0, (10) 

with initial condition M = I, the NxN identity matrix. Although this equation 
is linear, it is coupled with Eq. (9) through the dependence of the Jacobian ma- 
trix F on the trajectory s t . Equation (10) is obtained by differentiating Eq. (9) 
with respect to s t{J and applying the chain rule. That there exists a unique 
solution M G R NxN of Eq. (10), in fact with one continuous time derivative, 
for as long as the trajectory s t exists in 5, is guaranteed by Assumption 2 and 
the linearity of Eq. (10). The dependence of the solution M on the trajectory 
s t can be expressed fully as M = M(t) = M(f;s to ), since s to determines the 
trajectory s t . 

Now let M tjto denote the unique solution of Eq. (10) in R iVxJV , over a finite 
time interval for which s t exists in S, that corresponds to the initial condition 
M t0jto = I. From the preceding discussion it follows that for each point q to G 
M. N , the linear equation 

^+F(s t ,f)q = 0 (11) 

has unique solution 

q t = M t)to q to G 1 W , (12) 
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and that q t has one continuous time derivative, for as long as s t exists in S. The 
linear equation (11) is called the (deterministic) perturbation equation associ- 
ated with the original nonlinear dynamics (9). The matrix M( j(o = M t! t 0 (s to ), 
which according to Eq. (12) expresses the solution of the perturbation equation 
directly in terms of its initial condition, is called the fundamental matrix, or so- 
lution operator, of the perturbation equation. The analogue of Eq. (11) for the 
stochastic initial value problem, which gives the evolution of stochastic initial 
perturbations under second-moment closure, is derived in Sec. 5 along with the 
corresponding mean and covariance evolution equations. Energetic consistency 
of the mean and covariance evolution equations is demonstrated in Sec. 6. The 
fundamental matrix of the stochastic perturbation equation has special proper- 
ties due to this energetic consistency, which are described in Sec. 9. 

It is well known (e.g. CODDINGTON and LEVINSON [1955, Ch. 1, Thm. 7.3]) 
that Eq. (10) can be solved explicitly for the determinant of M t to : 


det M Mo = exp 


f trF(s t ,t)cLt 
Jt . 0 


(13) 


where tr F denotes the trace, or sum of the diagonal elements, of F. Define the 
symmetric and anti-symmetric (skew-symmetric) parts of F as F s = i(F + F r ) 
and F a = i(F — F^ ), respectively, where the superscript T denotes transpo- 
sition, so that F = F s + F a and F 1 = F s — F“. Since the diagonal elements 
of a real skew- symmetric matrix are all zero, F can be replaced in Eq. (13) by 
the symmetric matrix F s . This is one indication of the important role that the 
symmetric part of the Jacobian matrix plays in the dependence of trajectories 
on their initial points. In Sec. 6 it will be seen that the exchange of energy be- 
tween the mean state and the stochastic perturbations occurs solely through F s . 
This role of F s in the energetic coupling of the mean and covariance evolution 
equations is examined in detail in Secs. 7 and 8. 

Another immediate consequence of Assumption 2 is that if di /dt £ C(SxT), 
then in fact s t has two continuous time derivatives, not just one. This follows 
by differentiating Eq. (9) once with respect to time and applying the chain rule: 


cP s 
dt 2 


F(s, f)f (s, t) 


<9f(s, t) 
dt 


Also, if Of/dt £ C(S x T), then Assumption 2 implies Assumption 1. It will 
not be assumed that di/dt £ C(S x T), although this often does hold. For 
instance, in many problems f does not depend explicitly on time. 


3 Conservation of Total Energy 

Suppose now that one is given a nonlinear system 


dx 

dt 


+ g(x) = 0, 


(14) 


7 



with state space X, and with g : X x T — > K w satisfying the hypotheses 
(Assumptions 1 and 2) imposed earlier on f . Suppose also that there is a known 
function s : X — > M. N such that the “total energy” 

E(t) = ±s T (x t )s(x t ) (15) 

is conserved by the solutions of Eq. (14), that is, 


dE_ 

dt 


s T (xt) 


dsjpct) 

dt 


= 0 , 


(16) 


for each trajectory x t . Suppose finally that s = s(x) defines a continuously dif- 
ferentiable coordinate transformation between the state space X and the range 
s(A), and denote by A(x) = ds/dx the Jacobian matrix of this transformation. 
Then for s = s(x t ) one has 


ds 

dt 


a( x ; 


dx 
dt ’ 


and so from Eq. (14) it follows that s(x t ) satisfies the nonlinear system (9), for 
f given by 

f (M) = A(x(s))g(x(s),t) , 

where x = x(s) is the inverse transformation of s = s(x). Furthermore, for 
Eq. (9), the total energy becomes simply E = Js T s, and the statement (16) 
that the total energy is conserved becomes simply s T f = 0. Thus, rather than 
considering Eq. (14) directly, for some general energy expression, in this article 
Eq. (9) is considered instead, under the simple hypothesis that the total energy 
E = ^s T s is conserved: 


Assumption 3 For all s G S and t G T, 

s T f(s, t ) = 0. 


(17) 


Transforming a general energy expression out of the problem in this way simpli- 
fies substantially the study of energetics of the second-moment closure equations. 
The cost is potentially a complicated expression for f. 

Assumption 3 says that the trajectories s t of Eq. (9) satisfy 

T T 

s t s t = s t o S *o> 

for as long as t G T and the trajectory exists in S. In geometrical terms, 
this means simply that each trajectory remains on the hyperspherical surface 
s T s = 2 E in on which it originates at time to- A trajectory can cease to 
exist only if this surface intersects the boundary of S. 

Since the change of coordinates to “energy variables” s is central to the re- 
sults of this article, it is worthwhile to consider how it works in simple examples. 
Consider first the quadratically nonlinear system 


— - + cu = 0, 
at 




with c a nonzero constant. Assumptions 1 and 2 are satisfied with the state 
space taken to be all of R 2 , and with arbitrarily large final time T for the 
interval T = [fo,T]. Therefore, for each initial condition in R 2 , this system 
has a unique solution in R 2 which exists over arbitrarily long time intervals. 
In addition, for E = ^(u 2 + <j>) one has dE/dt = 0. This “energy” expression 
suggests that one consider also the state space A consisting of the upper half- 
plane <f> > 0 in (u, 0)-space, so that E > 0. The change of variables Si = u, 
S 2 = 4 )1 / 2 is a continuously differentiable coordinate transformation from X 
onto itself, and it yields Eq. (9) with 



which is singular along the Si-axis S 2 = 0. Assumptions 1 3 are satisfied by this 
f, with S being the upper half-plane S2 > 0 in (si, S2)-space, and again with 
arbitrarily large final time T for the interval T. The “total energy” E = is T s is 
conserved on each trajectory of Eq. (9), and every trajectory exists either until 
it hits the Si-axis, which is the boundary of S , or for all t £ T if it never hits 
the si-axis. 

The solution of Eq. (9) in this example, for each s to £ S , is just 
si(t) = Si(f 0 )e -c(t-to) , 
s 2 (t) = [2 E - sl(t)] 1/2 . 

There is one stationary solution, which is Si = 0, S 2 = ( 2E ) 1 / 2 . If c > 0 
and si (to) 7^ 0, then |si(t)| decreases monotonically toward zero, and S 2 (t) 
increases monotonically toward (2E) 1 / 2 : all solutions approach the stationary 
one, and no trajectory can ever hit the Si-axis. Thus if c > 0, then the trajectory 
corresponding to arbitrary s to £ S exists in S over arbitrarily long time intervals. 
If c < 0 on the other hand, and if Si(<o) 7^ 0, then |si(t)| increases monotonically 
toward (2 E) x ! 2 and S2(t) decreases monotonically toward zero, and the solution 
ceases to exist at time T*, 



when s 2 (T*) = 2 E, S 2 (T*) = 0, and f is no longer defined. However, for the 
problem in the original variables u and 4>i as posed on all of R 2 , nothing bad 
can happen at this time T*, nor at any other finite time, since it was already 
seen that every trajectory of the original problem exists in R 2 over arbitrarily 
long time intervals. In fact, all that happens is that the trajectory continues 
downward into the lower half-plane 4> < 0 in (u, 4>)~ space, along the parabola 
4>t = 2 E — it 2 . Thus it is possible for long-time existence of solutions to be lost 
through the transformation to energy variables. 
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Fortunately, for spatially discretized versions of the main first-order hyper- 
bolic partial differential equations of atmospheric and ocean dynamics, nothing 
need be lost in the transformation to energy variables, since all that is required 
in the transformation is to take the square root of mass (gravitational poten- 
tial energy) and temperature (internal energy) variables that are required on 
physical grounds to remain positive. As a simple example, consider the shallow- 
water equations in one space dimension, taken to be periodic. These are usually 
written as the momentum equation 


du 

dt 


du dcj) 
U dx dx 


0, 


and the continuity equation 

(ty dufi 
dt dx 

where u is the speed, (j> is the geopotential, and x € [0, L\ is the space vari- 
able. From the differential equations it follows that solutions satisfy the energy 
equation 

d(4>u 2 + <j) 2 ) d(cj)u 3 + 2 (j) 2 u) 

Oi + d~x = 

which on the periodic domain [0, L\ implies conservation of the total energy 


E = 



(4>u 2 + </> 2 ) dx. 


Reasonable spatial discretization of the variables u and <fi gives rise to a system 
of ordinary differential equations (14) that conserves a discretized version (15) 
of this energy integral. 

The change of variable from u to a = <ft~! 2 u transforms the momentum 
equation to 


da da du 

~m +u ^ + ^di a 


+ 0 1/2 ^ = o, 

OX 


and substituting u = </> -1 / 2 a here and in the continuity equation gives the 
shallow- water system in terms of the variables a and <f> alone. Reasonable 
spatial discretization of the transformed system, say with a n (t) = a(x n ,t), 
4>n{t) — <f>(x n ,t), and x n = nL/N, for n = 1, ...,1V, gives a system (9) 
for s = [«i, . . . , q:jv ; 0i, • • • , 0 jv] T that conserves the discretized total energy 
E = ^s T s-^-. One can see that, because of the substitution u = 0 _1 / 2 a, the 
function f will be singular on each of the hyperplanes (j> n = 0, n = 1, . . . , N, in 
R 2W . Merely taking S to be the convex open set 


<5 = {s = [ai, . . . , «at, (j) i, • • ■ , <I>n] ■ (frn > 0, n = 1, . . . , N} 


does not, however, give rise to an initial-value problem whose trajectories exist 
in this S over arbitrarily long time intervals. That is, trajectories starting in 
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this S may hit the boundary of <S, developing zero geopotential somewhere, at 
some time T* < T, and thereby cease to exist. 

For the one-dimensional shallow-water equations, it is straightforward to 
define a state space such that every trajectory does exist over arbitrarily long 
time intervals, for many spatial discretizations of either the (u, (f>) system or the 
(a, 4>) system. It follows from the shallow-water equations that the characteristic 
speeds 

c + =u + 0 1 / 2 = </> -1 / 2 (a + (j > ), 
c_ = u — q i 1 / 2 = </> -1 / 2 (a — 0), 
satisfy the coupled advection equations 


dc+ ,a t x dc + 

— + ( 3 C + + 3 C_)— - 0, 


dc _ 
dt 


dc _ 


+ (\c + + ±c-)'- dx 


= 0. 


Therefore if initially c_|_ > 0 and c_ < 0, equivalently <j> > |a|, for all x G 
[0, L\, then this remains true for all time. Positivity-preserving (of c+ and — c_) 
discretizations (e.g. LIN et al. [1994]) of the shallow- water equations therefore 
have the property that trajectories with initial points in the convex open set 


Sq = {s = [ai, ... ,cxn, (/>!,-■■ ,(/>n] ■ (t>n > 


= !,•••, N} 


exist in So over arbitrarily long time intervals. More generally, if there is a con- 
stant c > 0 such that initially c+ > c and c_ < — c (equivalently > c+ |it|, 
or t/) 1 / 2 > c/2+(c 2 /4+|a|) 1 / 2 ) for all x € [0, L], the coupled advection equations 
imply that this also remains true for all time. Therefore, for appropriate spatial 
discretizations, trajectories with initial points in the convex open set 

S c = {s = [a\, . . . , aiy, <j> i, . . . , 4>n] ■ 4\! 2 > c/2+(c 2 / A+lcxnl) 1 ^ , n = 1, . . . , IV} 

exist in S c over arbitrarily long time intervals. On the state space <So; the 
geopotential can approach zero arbitrarily closely wherever u = 0. On the state 
space S c with c > 0, however, the geopotential is guaranteed to be strictly 
positive everywhere, (f) > (c + |w|) 2 > c 2 , for all t G T, and therefore the total 
energy is also strictly positive, 

nL pL 

E = \ / (a 2 + </> 2 ) dx > | / q i 2 dx > \Lc A . (18) 

Jo Jo 


4 The Stochastic Initial- Value Problem 

Now let be the sample space of a complete probability space, and denote by £ 
the expectation operator on the probability space. A random vector s G is a 
vector function s = s(u>) of the probability variable u> £ 12, that is, s : 14 — * i w . 
For fixed u>, s (u>) is called a realization of the random vector s. 
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Consider a second-order random vector, that is, a random vector s £ R w 
such that £s T s < oo. By the Schwarz inequality one has 

(£s) T (£s) < (£|s|) t (£ |s|) < £s T s, 

and therefore the mean s = £s exists and is finite. Denote the departure from 
the mean by s' = s — s. Since £s' = 0 one has 

£s t s = s T s + £s ,t s', 

and therefore the total variance V defined by 

V = £s ,t s' 

is also finite, V < £s T s. Since the total variance is finite, it follows by a further 
application of the Schwarz inequality that the elements of the N x N covariance 
matrix P defined by 

P = £s's’ t 

also exist and are finite. The covariance matrix is, by definition, symmetric and 
positive semidefinite. Also, the total variance is just the trace of the covariance 
matrix, V = trP. Associating randomness with uncertainty, one can say that 
the total variance of a second-order random vector is a scalar measure of its 
uncertainty. 

For a random vector s £ R w , one writes s £ S wpl (with probability one) if 
s (u>) £ S for all u> £ O, except possibly for an w set of probability measure zero. 
A random vector s £ S wpl will be called a random state. Every second-order 
random state s has mean s £ S, because S was taken to be convex. 

In the stochastic initial-value problem treated here, Eq. (9) is considered for 
second-order random initial states s to , that is, random initial states with 

£E to = ±£sf o s to < oo. 

Denote by s t (w) the trajectory corresponding to a realization s to (cu) £ S. From 
Assumption 3 it follows that s t (oj) satisfies the total energy conservation equa- 
tion 

sf(w)s tM = s[ 0 (u)s to (uj), (19) 

for each w £ fl such that s to (w) £ S. If the deterministic initial- value problem 
is such that all trajectories starting in S exist in S over arbitrarily long time 
intervals T = [to,T], for instance if S = R w , then the realizations s t (u>) with 
s to (w) € S define a random vector s t £ S wpl, for all t £ T, and by taking 
expectations in Eq. (19) it follows that s t is also a second-order random vector, 
for all t £ T . In this case, therefore, s t has finite mean s t £ S and covariance 
matrix P f £ R ArxJV , for all t £ T. The family of random vectors {s t : t £ T} 
is called a (second-order) stochastic process. The trajectories s t (w), for fixed 
u> £ £l, are called the sample functions, or sample paths, of the process. By 
Assumption 1, each sample path with s to (w) £ S is in C l (T). In other words, 
the sample paths are in C 1 (T) wpl. 
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In case the deterministic initial-value problem is such that there are trajec- 
tories starting in S that hit the boundary of S before time T, one can restrict 
the problem to initial states in some prescribed subset S' C <S, not necessarily 
convex or open, but whose boundary nowhere touches the boundary of S. Since 
the trajectories are continuous paths in state space, traversed at finite speed, 
there is then a time T' > to such that the trajectory s t (u>) corresponding to 
an arbitrary point s to (w) G S' exists in S for all time t in the closed interval 
T' = [t 0 ,T'], with T' independent of s to (w). Thus in this case, the stochastic 
initial- value problem is restricted to second-order random initial vectors s to G S' 
wpl. The realizations s t (w) with s to (w) € S' still satisfy Eq. (19), but now only 
for t G T'. Thus they define a second-order stochastic process {s t : t G T'}, with 
s t £ S wpl, and with finite mean s t G S and covariance matrix P t G R JVxJV , 
for all t £ T' . The sample paths are in C l (T') wpl. 

Thus two cases have been distinguished for the stochastic initial- value prob- 
lem. In the first one, Eq. (9) defines a second-order stochastic process {s t : i 6 
T} for each second-order random initial state. In the second case, Eq. (9) defines 
a second-order stochastic process {s t : t G T'} for each second-order random 
initial vector s to such that s to G S' wpl. Denote by V) the total variance 

V t = Ss' T s' t = trPf, 

for t G T in the first case, and for f G T' in the second. Taking expectations in 
Eq. (19) gives 

sj s t + V t = sJ o s to + V t0 , (20) 

for t G T in the first case, and for t G T' in the second. The total variance 
V t is a scalar measure of the uncertainty present in the random state s t due to 
uncertainty in the random initial state s to . 

Equation (20) says that the uncertainty in solutions of Eq. (9) due to uncer- 
tainty in the initial condition, as measured by the total variance, can increase 
(decrease) only as a result of extracting energy from (inserting energy into) the 
mean state s, with the change in total variance balanced exactly by twice the 
change in total energy ^s T s of the mean state. This is purely a consequence of 
conservation of total energy for the nonlinear dynamics, and it holds regardless 
of any assumptions one might make on the form of the probability distribution 
of the random initial state s to , apart from existence of its first two moments. 
In fact, Eq. (20) holds even if no moments beyond the first two exist at any 
time. It is simply a statement about second-order stochastic processes whose 
realizations satisfy Eq. (19) with probability one. 

Equation (20) implies in particular the bound 

V<sf 0 s t0 + V t0 , (21) 

valid for all time t G T in the first case, t G T' in the second, with equality 
holding at some particular time r if, and only if, all the energy of the mean state 
has been extracted at that time, s r = 0 . This is a simple, general statement 
of the maximum level of uncertainty that can occur in solutions of Eq. (9). 
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There is no implication, however, that the uncertainty actually saturates, that 
is, becomes or approaches a constant in time, either at the level sJ o s to + V to or 
at any other level. Saturation of uncertainty will be discussed briefly at the end 
of Sec. 8. 

For many dynamical systems, the state space is naturally bounded away 
from the origin of coordinates s = 0 in R iV , so that in fact the mean state 
cannot vanish at any time and the bound (21) can be improved upon. This 
is the case when there are mass-like and/or temperature-like state variables, 
for instance, that are constrained by the physics of the problem to be bounded 
from below by positive constants. If every point s in the state space S of the 
problem satisfies an inequality s T s > 2 E m - ln > 0, then since s t £ S one has 
sjs t > 2E m - m > 0; here 2£’ m ; n is just the minimum Euclidean distance from the 
origin to the boundary of S. In this case Eq. (20) implies the stronger bound 

Vt < s to S *o + Vt 0 — 2 .Emin- (22) 

For example, in Sec. 3 it was shown that for appropriately discretized versions 
of the one-dimensional shallow-water equations with state space S c , c > 0, the 
geopotential satisfies > c 2 for all t £ T, and since the discretized total energy 
was E = ^s T s-^-, Eq. (18) gives sjs t > Nc 4 . Thus for appropriate discrete 
shallow- water dynamics on S c one has 

V t < sf 0 s to + V to - Ac 4 , (23) 


for all t £ T. 

Since Eq. (20) holds without any assumptions on moments beyond the first 
two, even without their existence, one expects to be able to find approximate 
evolution equations for just the mean and covariance matrix, such that their 
solutions are guaranteed to satisfy Eq. (20) and therefore also the bounds (21) 
and (22). A closed system of differential equations for the mean and covariance 
matrix that has this property will be said to be energetically consistent. The 
second-moment closure equations derived in Sec. 5 constitute a closed, non- 
linearly coupled system of differential equations for the mean and covariance 
matrix. In Sec. 6 it will be shown that the nonlinear coupling term in the 
second-moment closure equation for the mean makes them energetically consis- 
tent. First the exact equations for the mean and covariance matrix, which are 
not closed unless the dynamics are linear, will be derived. 

If E represents a physical total energy, then for both the deterministic and 
stochastic initial-value problems there seems little reason to consider initial 
points in R w with total energy greater than or equal to some prescribed, fixed 
amount, say E > E m ax . Such points will be eliminated from consideration by 
making the following simple hypothesis on S: 

Assumption 4 S C <S max , where <S max denotes the interior of the hypersphere 
s T s = 2 £ max in R^: 

S max = {s£R N : s T s < 2 £ max }. 


14 



This assumption imposes no restriction on existence of solutions, since no tra- 
jectory starting in S C S max can ever hit the boundary of iS max at any time, by 
conservation of total energy. It is, however, a restriction on the initial probability 
distribution, and therefore on the probability distribution at any time, because 
it implies that if s is a random state, then s T s < 2 E max wpl. In particular, 
no random state is (multivariate) normally distributed, and the marginal distri- 
butions of a random state also cannot be normal. This may seem a significant 
restriction. However, it is suggested by the physical problem at hand. Moreover, 
as discussed already, for atmospheric and ocean dynamics there are also usually 
state variables that are constrained to be positive, in fact often bounded from 
below by positive constants, and these cannot be normally distributed. 

An immediate consequence of Assumption 4 is that every random state s is a 
second-order random vector, in fact, £ s T s < 2 E max . Also, since f € C(SxT) by 
Assumption 1, it follows from Assumption 4 that f T f < oo on S x T . Therefore, 
if s is a random state, then f(s) is a second-order random vector, £f T (s)f(s) < 


Again let s t (either for t € T or t G T' . depending on the case) denote the 
solution of the stochastic initial- value problem for Eq. (9). Then £|f(s t , t)| < oo 
since f(s t ,t) is a second-order random vector, and therefore 


[ f|f(s T) 

Jt 0 


r) I dr < oo. 


It follows (e.g. DOOB [1953, Thm. 2.7, p. 62]) that 


£ f f(s r ,r)dr 
Jto 


£ f (s T 


I dr. 


(24) 


where both integrals exist and are finite. Now write Eq. (9) as the integral 
equation 

s t (w)=s toM- f f(s r (w), r) dr. (25) 

d to 

Taking expectations and using Eq. (24) gives 


s* — S( 0 


f(s T , r ) dr , 


'to 


(26) 


where f(s r ,r) = ff(s r ,T), which shows that the mean state s t is a continuous 
function of time. Differentiating in Eq. (26) gives an equation for the mean 
state in differential form, 

Jo 

■^+f(8 t) i) = 0, (27) 

which is satisfied almost everywhere in T in the first case, that is, except possibly 
on a subset of T of Lebesgue measure zero, and almost everywhere in T' in the 
second. This exact equation for the mean state is not a differential equation 
unless f(s,f) is linear in s, in which case = f(s t ,f). Also, f(s t ,f) is 
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not necessarily continuous in time, and so s t is not necessarily continuously 
differentiable. Comparing Eqs. (9) and (27) shows that the commutation 


c ds t _ ds * 
fit fii 


(28) 


holds, almost everywhere in time. 

To derive the equation for the covariance matrix, first subtract Eq. (26) from 
Eq. (25), to obtain 


= s t 0 (^) 


f'(s r (w),r)fir, 


'to 


(29) 


where s' = s — s and f' = f — f . Postmultiplying this equation by the transpose 
of itself, then taking expectations, and then exchanging the order of expectation 
and integration, which again can be justified by Assumption 4, gives 

Pf = Pfo - [ £[f'(Sr,T)s Z]dT-[ £[s' t J ,T (s Tl T)} dr 

Jto Jto 

+ f f £ [f'(s ri ,Ti)f' T (s r 2 ,T 2 )] dn dr 2 , 

Jto Jto 

which shows that the covariance matrix P t is a continuous function of time. 
Differentiating this result and using Eq. (29) gives an equation for the covariance 
matrix in differential form, 

^ + £ [f'(s t , t ) sf ] + £ [f'(s t , t)s[ T ] T = 0, (30) 

satisfied almost everywhere in T, or in T . Again, this is not a differential equa- 
tion unless f is linear, and the elements of Pf are not necessarily continuously 
differentiable. Equation (30) can be used to show that the commutation 


c ds t s t T = 

dt dt 


(31) 


holds, almost everywhere in time. 

Equations (27) and (30) are usually derived under an hypothesis much 
weaker than Assumption 4 which, like Assumption 4, also implies that f(s t ,f) 
is a second-order random vector (e.g. DOOB [1953, p. 277, hypothesis H 2 ], 
JAZWINSKI [1970, p. 105, hypothesis Hi]). However, Assumption 4 makes 
sense for conservative dynamics, and it greatly simplifies the derivation of Eqs. (27) 
and (30). It also makes for little difference between the formulation of the 
stochastic initial-value problem and that of the deterministic initial- value prob- 
lem, since it makes every random state a second-order random vector. All that 
has been necessary for the stochastic problem was to restrict initial random 
vectors to lie in some set S ' (wpl) contained wholly in the interior of S in case 
solutions of the deterministic initial-value problem do not exist over arbitrarily 
long time intervals, to ensure that all the sample paths exist for some minimum 
amount of time T' — to > 0 wpl. 
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5 The Second-Moment Closure Equations 

Two final hypotheses on f will now be made: 

Assumption 5 <9 2 f(s ,t)/dsjdsk € C(S x T), for j,k = 1 


Assumption 6 The second partial derivatives of f are Lipschitz continuous in 
s on S x T, uniformly in t. That is, there are constants Kjk such that 


d 2 f(si ,t) 
dsjdsk 


<9 2 f(s 2 ,t) 
dsjdsk 


< Kjk ||si 


S 2|| , 


where || ■ || denotes the Euclidean norm on M. N , for each si,S 2 £ S, for each 
t £ T, and for j, k = 1, . . . , N. 


Let s be a random state, that is, s £ S wpl. Since s = £s £ S, it follows 
from Assumption 5 that f(s) = f(s, t) has a Taylor expansion about s up to 
second order, 


f(s) = f(s) + F(s)s' + \ ds^dsk 3 '^^ + ^ ( 32 ) 

j k 3 

where F = F(s) = F(s ,t) is the Jacobian matrix of f introduced in Sec. 2, 
and where s' is the j th element of the random vector s' = s — s. It follows 
from Assumption 6 that the remainder term g is 0(( s') 3 ) for each fixed s £ S. 
It followed from Assumption 4 that s is a second-order random vector, and 
along with Assumption 1 that f(s) is also a second-order random vector, so that 
f = £{ exists and is finite. Taking expectations in the Taylor expansion then 
shows that g = £g also exists and is finite, and that 

f = f(s) + i^^^Jp ife +g, (33) 

j k 3 


where Pjk = £ s'js' k is the (j, k) th element of the covariance matrix P = £s's' T . 

The mean equation is obtained by substituting Eq. (33) into Eq. (27) and 
neglecting the remainder term g, to yield 


ds 

dt 


+ f ( s ) + 

j k 


d 2 m „ 

dsjdsk 


= 0 , 


(34) 


where the time subscripts have been omitted because the mean equation is 
only an approximate equation for the evolution of s t . The covariance evolution 
equation is obtained by using Eqs. (32) and (33) to approximate f' = f'(s) = 
f(s) — f by F(s)s', and substituting this into Eq. (30), yielding 

dP 

— +F(s)P + PF t (s) = 0, (35) 
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where the time subscripts have again been omitted because the covariance evo- 
lution equation is only an approximate equation for the evolution of P 4 . 

Equations (34) and (35) together constitute the second-moment closure equa- 
tions for the evolution of the mean state and covariance matrix. In case f = 
f(s ,t) is linear in its first argument, the second-moment closure equations de- 
couple and they are exact: the Jacobian matrix F = F(s, t) is independent of 
its first argument, so that the covariance evolution equation decouples from the 
mean equation, and the second partial derivatives of f (s, t) with respect to their 
first argument all vanish, so that the mean equation likewise decouples from the 
covariance evolution equation. In case f is nonlinear, the mean and covariance 
evolution equations are fully coupled, in both directions. 

Let V C R NxN be an open set, and suppose that s to G S and P to G V. Then 
there is a half-open time interval Tq = [to, To), with To depending in general 
on s to and P to , such that there exists a unique solution (s (t) € <S, P(i) G V) of 
this coupled system for all t G %, satisfying initial condition (s(t 0 ), P(to)) = 
(s to ,P to ). Assumptions 1, 2, and 5 guarantee that there exists at least one solu- 
tion on 7o satisfying the initial condition, and also that every solution on 7o is in 
C ,1 (7o). Assumptions 2 and 6 guarantee that there exists at most one solution 
on 7o satisfying the initial condition. As in Sec. 4, if s io and P to are restricted to 
lie in some subsets S C S and V C V, respectively, whose boundaries nowhere 
touch the boundaries of S and V, then existence and uniqueness of solutions 
are assured over some closed time interval T = {to, T] with T > to. Thus if one 
is interested in solving Eqs. (34) and (35) with some particular s to G S given, 
as is often the case, then existence and uniqueness on some closed time interval 
T = [to,T] are assured for any P«o e V, simply by taking S = {s to }. 

In fact, Eqs. (34) and (35) are supposed to be solved for initial condition 
(s to ,P to ) being the mean state and covariance matrix of a random state s to , 
if the equations are to approximate the evolution of the mean and covariance 
matrix of the random state s t . In particular, P <0 is supposed to be symmet- 
ric positive semidehnite. Moreover, this means that P to cannot be specified 
independently of s to . For instance, under Assumption 4 one has the simple 
restriction 

sJ o s to + tr P t 0 < 2 E max . (36) 

The particular geometry of S for a given problem can restrict P to in terms of 
s to much more severely than this, with the restriction being the strongest when 
s to is near the boundary of S, as discussed further below. Since Pt 0 is not 
supposed to be specified independently of s to , it is simplest to pose the initial- 
value problem for a given, fixed s to G S, then to define the set V = V(s to ) of 
symmetric positive semidehnite matrices over which P io is allowed to vary due to 
the particular geometry of S, and finally to define a subset V = V(s to ) C V{s to ) 
whose boundary nowhere touches that of V(s to ). In this way, one has existence 
and uniqueness on a closed time interval T = [foiT], for the given s to G S 
and for all P to G V(s to ), along with the assurance that every P to G V(s to ) 
is the covariance matrix of a random state s (o with the given mean s to G S. 
Then one can conclude that a physically meaningful problem has been posed. 
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If one is given both s to £ S and any particular p i 0 e V(s to ), then existence 
and uniqueness are guaranteed on a closed time interval T = [to,T] by taking 

r 

As a simple example of how the requirement that P to is the covariance matrix 
of a random state with mean s to makes Pf 0 depend on s to , consider again the 
shallow-water system described at the end of Sec. 3. There it was shown that 
on the state space S c , solutions exist over arbitrarily long time intervals, for 
appropriate spatial discretizations. In terms of the variable u rather than a, the 
space S c is described by the inequality 

ul + 2c\u n \ + c 2 < <t> n , 

for n = 1, . . . ,N. Taking expectations gives 

£ {u' n f + u 2 + 2c£\u n \ + c 2 < </>„, 

and since \u n \ < £\u n \, this implies that 

£(u'n) 2 <fin-(\u n \+c) 2 . (37) 

Since the mean state is in S c , one has (|u n | + c) 2 < </> n , so that the right-hand 
side of inequality (37) is indeed positive. This inequality is a restriction on the 
variance £ (u' n ) in terms of \u n \ and (j) n for every random state on S c , and the 
restriction becomes stronger as the mean state approaches the boundary of 5 C , 
that is, as the right-hand side of inequality (37) becomes small. For given <p n , 
n = 1, . . . , N, the restriction is mildest when the mean state is a state of rest, 
u n =0 for n = 

Returning to the general problem, suppose now that the set V = V(s to ) 
restricted by positive semidefiniteness and the geometry of S has been defined 
for each s to £ S. From inequality (36) and the fact that S is an open set it 
follows that V(s t , 0 ) is an open set, for each s to £ S. A simple choice for the set 
V = V(s to ) C V(s to ) of initial covariance matrices P to is the set 

V = = {P G V ■■ ip G V}, (38) 

for any /.t with 0 < /i < 1, which will be made for the sake of definiteness in 
Sec. 9, where a specific choice of V is needed. Since V(s to ) is an open set, 
Vfj(s to ) is also an open set. Note that both V{s to ) and ^(s t „) contain the 
origin in R NxN , for every s to £ S, since S is an open set. 

In case the solution of the deterministic initial-value problem for Eq. (9) 
exists over arbitrarily long time intervals, one would like to know whether or 
not the solution of the initial-value problem for Eqs. (34) and (35) also exists 
over arbitrarily long time intervals, for each s to £ S and P*o e P(s t0 ). This 
important question is not addressed in the present article. Also, the question 
of how well the solution of Eqs. (34) and (35) approximates the mean and 
covariance matrix of the stochastic process defined in Sec. 4 is not considered 
here. 
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For given s t(J G S and all p *o e V(s to ), there exists on some closed time 
interval T = [to, T] a unique solution of Eqs. (34) and (35) satisfying the initial 
conditions. From this point onwards, s t and P t (defined originally in Sec. 4) 
are redefined as this solution. The random vector s( = s [(to) is also redefined, 
for each to € fi such that s to (w) = s to + s( o (to) G S, as the unique solution on T 
of the (stochastic) perturbation equation 


ds 1 , 

s + F(f)s = °' 


(39) 


corresponding to given initial condition s( o (w) = s to (to ) — s to . That there exists 
a unique solution on T for each such to is guaranteed by Assumption 2, since 
the perturbation equation is linear. Finally the stochastic process {s t : t G T} 
is defined, as the one whose sample paths are given by s t (to) = s t + s((w). The 
mean state of this second-order process is s t , and P t is its covariance matrix. 
Note that for this process, 


ds t ds t ds( 

dt dt, dt 


(40) 


is continuous on T wpl, since ds t /dt is continuous on T and ds(/dt is continuous 
on T wpl. 

The stochastic perturbation equation has the form of the deterministic per- 
turbation equation (11) for the original nonlinear dynamics. When f is non- 
linear, the evolution of s( according to the stochastic perturbation equation 
depends through the Jacobian matrix on the mean state, and therefore also 
on the covariance matrix, since Eqs. (34) and (35) are fully coupled. Thus in 
the second-moment closure framework, to simulate the evolution of individual 
sample paths s t (to) requires first solving for not only the mean state, but for 
the covariance matrix as well. This is one consequence of the presence of the 
nonlinear coupling term in the mean equation. 

To close this section, it will be established that Eqs. (28) and (31) hold for 
the newly defined stochastic process. First, because the stochastic perturbation 
equation is linear, it has a fundamental matrix M ti(o , not to be confused with 
that of the deterministic perturbation equation defined by Eq. (10). The fun- 
damental matrix of the stochastic perturbation equation is defined for all t G T 
as the solution of the deterministic, linear equation 


dM t 't 0 

dt 


+ F( s )M*,t 0 — 0 


(41) 


corresponding to initial condition M t0jto = I, the N x N identity matrix. There- 
fore it expresses the solution of the stochastic perturbation equation directly in 
terms of random initial condition s( o : 

s* = Mt,tX 0 - (42) 


This fundamental matrix depends in general on the mean state, and therefore 
also on the covariance matrix, due to the nonlinear coupling term in the mean 
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equation. This dependence can be expressed fully as Mf jto = M tjto (s to , P to ), 
since s to and P to determine s t and P t uniquely, for all t G T . The fundamental 
matrix does not depend on the probability variable u>. Taking expectations in 
Eq. (42) therefore gives £s' t = 0 since £s' to = 0 , and then taking expectations 
in Eq. (39) gives 


since the Jacobian matrix does not depend on to. Taking expectations in Eq. (40) 
then yields Eq. (28). 

It can be verified that the same fundamental matrix can be used to express 
the solution of the covariance evolution equation (35), as 

P t = M t)t 0 P t 0 M£ t0 , (43) 


which is the operator form of the covariance evolution equation. Since P to = 
£s' to s^ is symmetric positive semidefinite, it follows from Eq. (43) that P* is 
also symmetric positive semidefinite. From Eq. (42) one has 


/ IT 


= M 


JT 


Mo s to s to ^Mo’ 


and on taking expectations it follows that 

Pt =£<*?■ (44) 

Postnrultiplying Eq. (39) by s' T , then adding the transpose of the result to itself, 
then taking expectations and using Eq. (44), leads to 

£ ^r + F ^ )P + pfT (^ = 0 

Comparing this result with Eq. (35) gives Eq. (31). 

To summarize: Attention will now be focused on the mean equation (34), 
the covariance evolution equation (35), which can be written equivalently in 
operator form as Eq. (43), and the stochastic perturbation equation (39), which 
can be written equivalently in operator form as Eq. (42). The covariance matrix 
can be taken to be defined by Eq. (44). The fundamental matrix is defined by 
Eq. (41). 


6 Energetic Consistency of the Second-Moment 
Closure Equations 

The second-moment closure equations were derived without reference to As- 
sumption 3 that total energy is conserved by the nonlinear dynamics (9). As- 
sumption 3 will now be used to show that the second-moment closure equations 
are energetically consistent. 
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By Assumption 5, Eq. (17) can be differentiated twice with respect to each 
of the N state variables and then evaluated for any s £ S and t £ T . Doing so 
once gives 

s T ^- + e fcT f = 0, (45) 

OSk 

for each s £ S and t £ T, and for k = 1, ...,1V, where e k denotes the k th 
column of the N x TV identity matrix. This can be written equivalently as 

f(s) = — F t (s)s, (46) 


which is a special relationship between f and its Jacobian matrix. Equation 
(46) implies in particular that if 0 £ S, then 


f(0) = 0, 


(47) 


which means simply that the nonlinear dynamics are not externally forced, 
and that s = 0 is a steady-state solution of the nonlinear dynamics. Recall, 
however, that typically 0^5 for geophysical dynamics, since there are usually 
state variables that must be positive on physical grounds. 

Differentiating Eq. (45) gives 


d 2f + e jT_^_ + e kT^_ 

dsjdsk dsk dsj 


(48) 


for each s £ S and t £ T, and for j,k = 1, ... ,N. The symmetric and anti- 
symmetric (skew-symmetric) parts of the Jacobian matrix, respectively F s and 
F a , were defined following Eq. (13). Equation (48) can be rewritten in terms of 
F s , as 


F* fe (s) 


it d 2f ( s ) 

2 dsjdsk ’ 


(49) 


for each s £ S and t £ T, and for j, k = 1, . . . , TV. Energetic consistency of the 
second-moment closure equations is an immediate consequence of this special 
relationship between the symmetric part of the Jacobian matrix and the second 
partial derivatives of f. 

To see this, premultiply the mean equation (34) by s T , and evaluate Eqs. (17) 
and (49) at s = s £ S, to find that 


x ds T s 
2 dt 


EE f M^ = °' 

j k 


(50) 


Recall that the total variance V is the trace of the covariance matrix, 

V = trP = £s' t s', 


where the second equality follows from Eq. (44). Applying the trace operator to 
the covariance evolution equation (35), and using the property that trPF T = 
tr F t P, gives 


1 dV 
2 dt 


+ tr F s (s)P = 0. 


(51) 
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This result follows also from the stochastic perturbation equation, as it must, by 
premultiplying Eq. (39) by s ,T , then using the fact that s' T (u)F“(s)s'(w) = 0 
for each to €E fi, since F a is skew-symmetric and s' is real, and then taking 
expectations. Adding Eqs. (50) and (51) gives 

= o, (52) 

at 

which is the statement (20) of energetic consistency for the second-moment 
closure equations. It implies in particular the bounds (21) and (22) on their 
solutions. 


7 The Role of the Symmetric Part of the Jaco- 
bian Matrix 

Consider for a moment the case of linear, conservative dynamics. If f(s,f) is 
linear in s, then the Jacobian matrix F = F(s ,t) is independent of s, and the 
second partial derivatives of f with respect to the state variables all vanish. 
Thus from Eq. (49) it follows that F s = 0, and therefore from Eq. (46) one 
has simply f(s) = F a s, with F a independent of s. The mean and covariance 
evolution equations, already simple for linear dynamics in general, simplify still 
further for conservative dynamics. 

That the Jacobian matrix has a symmetric part F s is one consequence of 
nonlinearity. Moreover, Eqs. (50) and (51) show that the exchange of energy 
between the mean state and the stochastic perturbations, which leads to the 
exact balance (52), occurs solely through the symmetric part of the Jacobian 
matrix. Equation (51) shows also that F s directly controls the growth and/or 
decay of the uncertainty measured by the total variance V. 

This section gives an overview of the role that F s plays in controlling the 
behavior of the total variance, and therefore in determining how the amount 
of energy exchanged with the mean state changes through time. Section 8 will 
then examine in more depth some of the special properties that F s has in this 
role, which are properties that result from conservation of total energy. 

According to Eq. (43), the solution P of the covariance evolution equation is 
symmetric positive semidehnite, with rank not exceeding that of P to . Therefore 
P has eigendecomposition 

P = WX 2 W T , 

where X 2 is the diagonal matrix of non-negative eigenvalues erf, . . . , er 2 , where 
rankP < L < N, and where W = [wj, . . . ,wj is the TV x L matrix of nor- 
malized real eigenvectors w i, wfw ; = 1 for l = 1 ,L. The eigenvalues and 
eigenvectors depend on time. It follows from the eigendecomposition that 

L 

F = trP = ^> 2 , (53) 

1 = 1 
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and that 


L 

trF s (s)P = ^a ; 2 w^(s)w ; . 
1=1 

Substituting into Eq. (51) gives 


E 


2=1 


2 dt 


w T F £ 


(s)w i 




= 0 . 


( 54 ) 


(55) 


Equation (55) says that the presence in P of an eigenvector w / with nonzero 
eigenvalue acts to increase (decrease) uncertainty if w^F s (s)w; < 0 (w^F s (s)w; > 
0). In particular, if F s (s,<) happens to be negative semidefinite, for all s £ S 
and t £ T, then the total variance V is monotone nondecreasing, dV/dt > 

0, and hence the total energy of the mean state is monotone nonincreasing, 
ds T s/dt < 0, independently of the initial condition (s to , P to ), for as long as the 
solution (s, P) exists. Thus dynamics with negative semidefinite F s constitute 
a worst case for predictability: the uncertainty V can only increase with time, 
or at best hold constant at times. That this should be a worst case is to be 
expected already from the deterministic initial- value problem: the deterministic 
perturbation equation (11) gives 

+ q T F s (sj,f)q = 0 

for perturbations of the deterministic trajectory s t , and so the size q T q of the 
perturbation is monotone nondecreasing, regardless of the initial perturbation, 
if F s is negative semidefinite on S x T . Similarly, if F s is positive semidefinite 
on S x T, then dV/dt < 0 and ds T s/dt > 0, independently of initial condition 
(s to ,P to ), for as long as the solution (s, P) exists. The dynamics in this case 
would be eminently predictable. 

By Assumptions 2 and 4, F s is bounded on SxT, and therefore the eigenval- 
ues of F s are bounded on SxT. The eigenvalues are real, since F s is symmetric. 
Denote the largest and smallest eigenvalues of F s by A max (F s ) and A m i n (F s ), 
respectively. In the next section it will be shown that, for conservative dynamics, 


A min (F s (s, f)) < 0 < A max (F s (s, t)), 


(56) 


at each point s £ S and t £ T. This means that F s cannot be either positive 
or negative definite, anywhere on SxT , although it does not rule out positive 
or negative semidefiniteness. Also it will be shown that at each point s in any 
open set 5* on which f is genuinely nonlinear at time t, as defined there, the 
inequalities in (56) are strict at time t: 

A min (F s (s, t))< 0 < A max (F s (s, t)), (57) 

which means that F s (s,t) has at least one positive and one negative eigenvalue 
at time t, at every point se5». Thus if F s is genuinely nonlinear on all of 5, for 


24 



all 1 6 T, then F s cannot be either positive or negative semidefinite, anywhere 
on S x T : there is potential for both growth and decay of uncertainty, at all 
times, regardless of where the mean state happens to be in S. Equation (55) 
shows in the genuinely nonlinear case that whether growth, decay, or neither 
actually occurs at a particular time depends on the eigenstructure of P, relative 
to that of F s (s), at that time. 

Now denote by u max = u max (F s (s, t)) and u min = u min (F s (s, <)), respec- 
tively, the eigenvectors of F s corresponding to eigenvalues A max = A max (F s (s, t)) 
and A m in = A m i n (F s (s, t)). It follows from Eqs. (53) and (54) that 

. trF s (s, f)P 

mm = Amin, 

P tr P 

where the minimization is over all symmetric positive semidefinite matrices P, 
and that the minimum is attained for P = P m ; n = 7 U m i n u^ in , with 7 an 
arbitrary positive constant. Furthermore, P m in € V(s) for all 7 small enough, 
where V is the set defined in Sec. 5, since V is an open set containing the origin 
in M. NxN . Also, P m in G P M (s) for each fi G (0,1), by taking 7 still smaller, 
where V ^ is the set defined in Eq. (38). Thus the minimum is achieved for 
matrices in 'P(s) and in V^s). Similarly, 

tr F s (s, t)P 

m&X 7 = ^maxj 

P tr P 

and the maximum is attained for P = P ma x = £>u max Um ax , with 6 an arbitrary 
positive constant. Rewriting Eq. (51) as 

,ldV trF«(s,*)P 
2 V dt trP 

then shows that — 2A m in (+2A ma x) is the maximum instantaneous relative rate 
of increase (decrease) of uncertainty, and is attained for P = P m in (P = Pmax)- 
The presence of both positive and negative eigenvalues of F s can lead to 
complex behavior for the total variance V as a function of time, and therefore 
also for the behavior of the total energy of the mean state as a function of 
time, which according to Eq. (52) mirrors precisely that of V/2. In general 
both growth and decay of total variance can occur, and therefore dV/dt can 
also vanish instantaneously. For instance, one can check that dV/dt = 0 at 

P — A m a X U m i n Umin A m in^max^ni ax * 

8 Genuine Nonlinearity and Essential Linearity 

Inequality (56) is readily established. First observe that Eqs. (17) and (46) 
together imply that 

s t F s (s)s = 0, (58) 
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for every s £ S and t € T, since F“ is skew-symmetric. Also note for later 
reference that if 0 £ <S, then Eq. (49) gives 

F s (0) = 0, (59) 

although again one should recall that typically 0 $ S for models of ocean and 
atmospheric dynamics. 

Equation (58) implies that at each point s 6 5, one has the following alter- 
native, at each fixed time t € T: either F s (s) has at least one positive and one 
negative eigenvalue, or else s is a null- vector of F s (s), 

F s (s)s = 0. (60) 


To see this, let 

F s (s) = U(s)A(s)U t (s) 

be the eigendecomposition of F s , where A is the diagonal matrix of eigenvalues 
Ai,...,Ajvj all of which are real since F s is real and symmetric, and where 
U = [ui, . . . , ujv] is the matrix of normalized real eigenvectors u;, iq = 1 for 
l = 1, . . . , N. Then Eq. (58) can be rewritten as 

N 

^Az(s)[s T iq(s)] 2 = 0. 

i=i 

Either all the terms in this sum vanish or there is at least one positive and one 
negative term, equivalently, at least one positive and one negative eigenvalue. 
The condition that all the terms vanish can be expressed as A(s)U T (s)s = 0, 
which is equivalent to Eq. (60) since U(s) is nonsingular. Thus the statement 
of alternatives has been demonstrated. 

Inequality (56) holds at each point s £ S, t € T where F s (s,t) has a null- 
vector. Therefore it holds at each s £ S, t, £ T where the second alternative 
condition, Eq. (60), holds. Furthermore, inequality (57) holds at each point 
s £ S, t £ T where the first alternative condition, that F s (s,t) has at least one 
positive and one negative eigenvalue, holds. Therefore inequality (56) holds for 
all s £ S and t £ T . The rest of this section will distinguish the two alternatives 
more tangibly, and examples will be given at the end of the section. 

The condition expressed by Eq. (60) can be expressed equivalently in a way 
that clarifies when it occurs and also allows one to check for its occurrence 
essentially by inspection of f, without even calculating F s . This is done by first 
introducing the polar coordinate p = (s 1 ^) 1 / 2 , so that for each s £ S one has 
s = pc with c on the unit hypersphere c T c = 1 in R w . Any scalar, vector, 
or matrix function <f) = (j>( s) = <j>(s, t) that is continuously differentiable with 
respect to s G S for all t £ T also has a continuous derivative with respect to 
p, for all s e S and t € T, and 

d(j) ^ d(j> dsi 
dp dsi dp 


V do 

^ d Sl Cu 
1 = 1 L 
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so that 


d (j> ^ d<p 

p lTp = ^~d7 l Sl ' 

By Assumption 2, </> can be taken to be f, in which case this relationship reads 

ms) 


p- 


d P 


= F(s)s, 


since <9f(s )/dsi is by definition the I th column of the Jacobian matrix F(s). 
Using Eq. (46), one then has 


F s (s)s= I [F(s)+F t (s)] s = \ 




= 2 P 


i-^ lf(3) (61) 


dp 


for every s £ S and t £ T. 

Thus the second alternative condition, Eq. (60), is equivalent to 


f(s) = p 


df(s) 

dp 


(62) 


It is always the second alternative that holds at the origin, if 0 £ S, as seen either 
by setting s = 0 in Eq. (60) or by setting p = 0 in Eq. (62) and using Eq. (47). 
Away from the origin, Eq. (61) says that the second alternative condition is 
equivalent to 


dp 1 f(s) 
dp 


= 0 . 


(63) 


Substituting Eq. (62) into Eq. (63) shows that, away from the origin, the second 
alternative condition is equivalent simply to 


9 2 f(s) 

dp 2 


= 0. 


(64) 


If Eq. (64) holds not just at a point s £ S, but for all points in an open 
set 5* C S, then f is linear as a function of p for all s £ 5,. By Assumption 
5, d 2 f/dp 2 is continuous on the state space S for each t £ T, and therefore if 
Eq. (64) holds on an open set <S* C S, then it holds also on 5* D S, where 5* 
denotes the closure of 5* in M iV . If f is linear in all the state variables, on an 
open set 5* C S, then certainly f is linear in p throughout <S*, and not only does 
the second alternative therefore hold on S *, but it was seen at the beginning 
of Sec. 7 that then in fact F s = 0 on 5,. It is possible for f to be linear in p 
on an open set 5, C5 without being linear in any of the state variables there. 
This is the case when f(s) = pg(c) for s £ 5*, for some function g whose first 
partial derivatives with respect to all the ci,l = 1, . . . , N, vanish nowhere for 
s = pc £ S*. 

In case condition (64) holds at every point s in an open set 5* = 5* (t) C S 
at a particular time t £ T, then f = f(s, t ) will be said to be essentially linear on 
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<S*(t). The equivalence between conditions (60) and (64), together with Eq. (46), 
implies that f(s , f) is essentially linear on an open set <S* (t) if, and only if, 

f(s, t) = F a (s, t)s, (65) 


for all s € <S*(f). The discussion at the beginning of Sec. 7 shows that if f(s, f) 
is linear in all the state variables on an open set <S*(t) C S, then Eq. (65) holds 
with F a (s, t) independent of s on <S*(t). Thus essential linearity, as defined here, 
amounts to generalizing the notion of linearity in such a way that one still has 
f = F a s, but with F a depending on s. 

In case condition (64) holds at no point s in an open set 5* = <S*(t) C S at a 
particular time t 6 T, then f = f (s, t) will be said to be genuinely nonlinear on 
5* (t) . With this definition, the original statement of alternatives can finally be 
rephrased, for open sets instead of points, as follows: F s (s,f) has at least one 
positive and one negative eigenvalue at each point s in an open set 5* (t) C S if, 
and only if, f(s,f) is genuinely nonlinear on <S*(t). It has also been shown that 
an open set 5* (t) C S on which f (s, t) is genuinely nonlinear cannot contain 
the origin s = 0, since Eq. (60) holds at the origin. If 0 ^ S, as is typical for 
geophysical problems, then it is possible for f to be genuinely nonlinear on all of 
S , for all time t € T. It is not difficult to show that for reasonable discretizations 
of the one-dimensional shallow-water system considered in Secs. 3, 4 and 5, f 
is genuinely nonlinear on all of the state space S c , for arbitrarily long time 
intervals T. 

To illustrate the ideas of genuine nonlinearity and essential linearity in the 
simplest setting, consider the case N = 2. It follows from Eq. (17) that f has 
the form 


f(s,f) = ft(s ,t) 


s 2 
-si 




for some scalar function ft, thus generalizing the first example of Sec. 3. It is 
immediate that, away from the origin, f is essentially linear precisely on those 
open sets <S*(f) on which /3 is independent of p at time t, that is, on which ft is 
a function only of the ratio S\/s 2 at time t. That Eq. (65) does indeed hold on 
every such set can be verified directly, by calculating F a . One finds that 


F“ 



1 

-1 0 ’ 


so that 


F“s= [ ft+hP^: 


dft 

dp 


s 2 

-Si 


and so Eq. (65) does hold where dft /dp = 0, as well as at the origin. Similarly, 
it is immediate that, away from the origin, f is genuinely nonlinear precisely 
on those open sets <S* (t) on which dft /dp vanishes nowhere at time t. The 
statement of alternatives says that on such a set, F s must have at least one 
positive and one negative eigenvalue, which for N = 2 means that F s must 
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have exactly one positive and one negative eigenvalue, and therefore must have 
a negative determinant, detF s < 0. One can calculate directly that 


det F s 


-Ip 2 



which verifies the statement. 

Finally, although it is not usually the case that 0 G S for geophysical prob- 
lems, it is instructive to consider dynamical behavior near the origin in case 
0 G S, particularly in light of the ideas of genuine nonlinearity and essential lin- 
earity. Recall from Eq. (47) that s = 0 is a steady-state solution of the original, 
deterministic nonlinear dynamics (9). 

Consider first the case in which f is essentially linear on an open set 5* (f ) C S 
with 0 G <S*(i). Thus f is linear in p on and hence f is linear in p near 

the origin along each coordinate axis. Therefore 

f k (ee l ,t) = a l k (t) |e|, (66) 


for some scalars a l k , k,l = 1 , ,N, and for all e small enough, where e l denotes 
the I th column of the N x N identity matrix. This follows from the fact that p = 
(s T s) 1//2 = |e| for s = ee l . By Assumption 2, dfk/dsi exists and is continuous 
at the origin, and since |e| is not differentiable at e = 0, it follows from Eq. (66) 
that a l k (t) = 0 for k,l = 1 Therefore f(s,t) = 0 for all s G <S*(t): 

essentially linear dynamics near the origin are trivial dynamics. Further, if f is 
essentially linear on all of S , for all t G T, and if 0 G S, what has just been 
shown is that then f = 0 on S x T. 

Now consider dynamics near the origin for general f. If 0 G S , one has that 
ee l G S, for l = 1, . . . , TV and for all e small enough. Evaluating Eq. (49) at 
s = ee l and using Eq. (59) gives 




lclT d 2 f(ee\t) 
2 dsjdsk 


and taking the limit here as e — > 0 then gives 


dF * fc (0,t) 

dsi 


1 d 2 fi(0,t) 

2 dsjdsk 


In view of Eq. (47), at s = 0 the mean equation (34) therefore reads 


dsj_ 

dt 


EE 


dF s jk (0,t) 

d Sl 


Pjk = 0, 


(67) 


for l = 1, . . . , TV. Thus the behavior of the mean state s at s = 0 depends on 
the symmetric matrices <9F S (0, t)/dsi, l = 1, . . . , N. In particular, even though 
s = 0 is a steady-state solution of the deterministic nonlinear dynamics, s = 0 is 
not necessarily a steady-state solution of the mean equation. The reason for this 
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is the presence of uncertainty in the stochastic dynamics, which is manifested 
through the covariance matrix in the nonlinear coupling term in Eq. (67). 

In view of Eq. (59), at s = 0 the covariance evolution equation (35) simplifies 
to 

dP 

— +F Q (0,t)P -PF“(0,t) =0. (68) 

Thus at s = 0 the covariance matrix evolves in an energetically neutral way, 
dV/dt = 0, and the evolution of the mean state at s = 0 still depends on 
this covariance evolution, through the nonlinear coupling term in Eq. (67). In 
particular, it is possible for the mean state to remain zero for a period of time, 
not just instantaneously, and subsequently to regain energy. That is to say, 
although it is possible for all of the energy of the mean state to have been 
extracted at some particular time t = r, and therefore for the upper bound in 
Eq. (21) actually to be attained at time r, it is also possible for the mean state 
to regain energy after time r, through interaction with the covariance matrix in 
the nonlinear coupling term. This effect is not present if the nonlinear coupling 
term is neglected. 

9 Bounds on the Growth of Relative Uncertainty 

It was shown in Sec. 7 that the instantaneous relative rate of change of the total 
variance V t satisfies the bounds 


— 2A max (F s (s i ,f)) < < -2A min (F s (s t ,f)), 

Vt at 

with each bound attainable by rank-one matrices P t G 7 ? M ( s t ) C V(s t ), for any 
given ft G (0, 1), where V and V are the sets defined in Sec. 5. It was shown 
in Sec. 8 that A m i n < 0 and A max > 0, and furthermore that if s t is in an open 
set on which f (s t , t) is genuinely nonlinear, then A m ; n < 0 and A max > 0. Both 
Amax and |A m i n | can be large, although finite, since no assumption has been 
introduced that would otherwise limit these values. Thus not only is it possible 
to have both growth and decay of total variance, at any instant of time it is also 
possible for the relative rate of growth or decay to be large. 

On the other hand, inequalities (21) and (22) show that V t cannot grow 
unboundedly. That is to say, although the total variance can grow rapidly at 
particular instants of time, growth over any interval of time is strictly limited, 
due to the energetic consistency of the mean and covariance evolution equa- 
tions which was demonstrated in Sec. 6. The present section examines the 
growth of Vt, relative to Vt 0 , over every interval of time for which the solution of 
the second-moment closure equations exists, thus providing time-independent 
bounds on the relative uncertainty Vj/V) 0 ■ The latter part of the discussion 
concerns inequalities (21) and (22) only, for the given state space S, without 
specific reference to the particular problem whose solutions satisfy the inequal- 
ities. Therefore the time-independent bounds obtained for V t /V to apply as well 
to the original stochastic process defined in Sec. 4, for as long as it exists. 
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Consider first the behavior of solutions of the second-moment closure equa- 
tions when the nonlinear coupling term in the mean equation is neglected, so 
that the fundamental matrix M tjto = M tjto (s to , P io ) introduced in Sec. 5 be- 
comes a function of s to alone, M f , to = M tjto (s to ). Denote by a t ,t 0 = cr tjto (s to ) 
the largest singular value of M t to , and denote by v (o = v to (s to ) the correspond- 
ing right singular vector, normalized so that v^v to = 1. The corresponding 
normalized left singular vector, u ( = Ui(s to ), is then given by 

u t = — M Mo v to . 

O'*, t 0 


Now take 

Pto = V to v to vf o , (69) 

so that V to = trP to , and use Eq. (43) to obtain 

Vt_ = trM t , to P tQ M^ o = tr a 2 tto V to \i t \iJ = 

V to tr P to tr P to 

The largest singular value <Jt,t 0 of any matrix M tjto also has the property that 


trM t|to PM^ o 

tr P 


< a 


2 

Mo » 


(70) 


for every symmetric positive semidefinite matrix P. This leads to the usual 
result that 

< <ri,t 0 (sto)> ( 71 ) 

v t 0 

with equality holding for the choice of P to given in Eq. (69). For this choice it 
is guaranteed that P to £ P M (s to ) C V(s to ), for any given /z £ (0, 1), by taking 
V to small enough. The bound in Eq. (71) depends on the initial mean state s to . 
Also, there is no general upper bound for the largest singular value cr* jto (s to ) 
itself. 

Now consider the behavior of solutions of the second-moment closure equa- 
tions with the nonlinear coupling term retained, so that the fundamental matrix 
is fully a function of P to , M tjto = M tjto (s to , P to ). The largest singular value 
of M t to is in general then a function of Pt 0 as well, cr t to = cr i to (s to , P to ), as 
are the corresponding right and left singular vectors. For any particular choice 
of P( 0 , for instance the one described in the preceding paragraph, inequality 
(70) still holds for every symmetric positive semidefinite matrix P. Therefore 
it holds for P = P to , and it follows that 


Vt_ 

Vt 0 


trM t , to P to M^ o 

tr P to 


— <T t,t 0 ( s *o>P*o)- 


Thus cr t 2 to is still an upper bound for Vt/Vt 0 , but the property that it is neces- 
sarily attained for some choice of P to has been lost. 
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What is actually desired is the least upper bound for V t /V to , call it of t (st 0 ): 

^2 ^ trM(, to (s to ,P (o )P to M^ t 0 (s t0 ,P t0 ) 

^Mo( s to) = sup — , (72) 

Pt 0 tr “to 

where the supremum is taken over all initial covariance matrices Pt 0 € P„K), 

for some given p G (0, 1), or over all P to in some chosen subset of P/i(st 0 ). It 

will be convenient to take the supremum over all initial covariance matrices in 
the open set P M ( s to ; V m in, V ma x) for some given fj,, V min and V max , where this set 
is defined for 0 < /i < 1 and 0 < V min < V max < 2 (£ max - E min ) by 

V^(s to ;V min , Vmax) = {P G P„(s to ) : V min < tr P < Vmax}- (73) 

This set is never empty, since P /t (s to ) is an open set. Taking the supremum over 
all P to G P/i(s to ; Vnin, Vmax) makes d t ,t 0 (st 0 ) defined in Eq. (72) depend also on 
Vnin and Vnax, that is, ( s^ 0 ) — uy£ 0 (s£ 0 , Vnin, Vnax)- Then one has 

VtfV t0 G; &t,to ( s *0 i ^niin, Vnax), 

with equality either holding, or holding arbitrarily closely, for some P to G 
P M (s to ; Vmin, Vnax)- One can also define 


® T (St 0 i Vnin , Vnax) — SUp &t,to (s^q , Vnin? Vnax), 
ter 


where T = [fo,T] C T is an interval of existence of solutions for all P to G 
Ef_i{st 0 ; Vnin, Vnax)- This yields 

V t /V t0 < a^(s t0 ; Vmi„ , V ma x) , (74) 


for all t G T, with equality either holding, or holding arbitrarily closely, at some 
time t G T and for some P to G P M (s to ; V m i n , V max )- 

The maximization problem of Eq. (72) is nonlinear, and in particular one 
cannot expect its solution to be independent of V to as it was seen to be in the 
linear case, when the nonlinear coupling term in the mean equation is neglected. 
It is for this reason that the supremum is to be taken over the set defined in 
Eq. (73), which allows a range V 0 G (Vnin, Vn a x) for the initial total variance. 
Inequalities (21) and (22) easily imply upper bounds for cr|^(s io ; V m in, V ma x), 
whose dependence on s to , Vnin and Vnax are explicit, as will be seen next. 

Consider first some cases in which 0 ^ S, so that inequality (22) holds. 
Rewriting it for V/V 0 gives 


Vt_ sJ o s to - 2E min 
Vt 0 V 0 


(75) 


Recall that s t(i and V 0 must also satisfy inequality (36), rewritten here as 


Vo < 2 E/max s t- 0 S to- 


(76) 
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The most unpredictable (largest possible V t /V to ) case occurs when s to is near the 
outer boundary of S, where the total energy is E max , for then Vt 0 in inequality 
(75) must be small, according to inequality (76). Therefore let sf o s to = 2 E max — 
e, with e > 0 fixed and small. Then inequality (76) holds if V to < e, so take the 
supremum over all Pio e V^(s to ',ae, e), for some given fi G (0, 1) and a G (0, 1). 
Setting sJ o s to = 2 E max — e and V) 0 > ae in inequality (75), and using inequality 
(74), then gives 


V t 2 

tt- < oWs to ;ae,e) <1 

Vt 0 a 


1 2 (E max — E m in) 


ae 


(77) 


for all t G T. This shows that Vt/Vt 0 can be large if e is small compared to 
2(A max — -Emin), even if a is taken to be close to one, that is, even if one allows 
only a small range V to G (ae, e) for the initial total variance. The upper bound 
(77) for the supremum ; ae, e) depends on the initial mean state s to only 

through its total energy E max — e/2. Recall from the derivation of inequality 
(22) in Sec. 4 that the bound (75) is tight if the mean state s ( corresponding to 
s to has energy near the minimum energy level E m i n . Therefore the bound (77) 
for <5^(s t 0 ', cte, e) is tight if some mean state with large initial energy E max — e/2 
approaches the minimum energy level E m in at any time t G T. 

The most predictable (smallest possible V t /V to ) case occurs when s to is near 
the inner boundary of S, where the total energy is E m j n , for then inequality (76) 
implies that V to need not be small. Let s^s to = 2E mm + e, with e > 0 fixed and 
small. Then inequality (76) holds if V to < Vmax = 2 (E max — E m , n ) — e, so take 
the supremum over all Pto € Vf,(s to -,aV max , Knax), for some given /u G (0,1) 
and a G (0, 1). Substitution into inequality (75) then gives 


A 

Vt 0 


< &y(s to ; aV m 


, km ax) < 


2(E max — E m i n ) + (k; — l)e 

2(£’ max — i?min) — e 


(78) 


for all t G T. This shows that Vt/Vt 0 can never be much larger than one if 
e is small compared to 2 (E max — unless a is also taken to be small, 

much smaller than e/2 (E max — E m i n ). The explanation for this is simply that 
V t < 2 (E max — E m i n ) for all t G T by Assumption 4, and so if V t(J is already close 
to the value 2(£' max — E m i n ), then there is no room for growth of the relative 
uncertainty Vt/Vt 0 . This behavior is far different than what occurs when the 
nonlinear coupling term is neglected, for then the initial total variance Vt 0 simply 
scales out of the problem. 

It is expected that in many applications E m i n and E max are both known 
reasonably well, with E max / E m i n not much larger than one. This still leaves 
plenty of room for growth of the relative uncertainty. Suppose the total energy 
of the initial mean state is the average of the minimum and maximum energy 
levels, so that s^s to = E m i n + E max . Then inequality (76) holds if V to < V max = 
Emax ~ Emin ■ Again taking the supremum over all P to G P M (s to ; aV max , Vmax), 
for some given p G (0, 1) and a G (0, 1), then gives simply 

—V < 4(s (o ;aV max , Vmax) < 1 + -, (79) 

Vt 0 a 
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for all t £ T. If the value of a is taken to be modest, say a = 1/3, then little 
growth of uncertainty can occur. However, if fo to is allowed to be small, say 
with a = 1/100, then considerable growth of uncertainty can occur, relative to 
this small value. 

In case 0 £ 5, the strict inequality (75) simply becomes non-strict, with 
Emi n = 0 also. The bounds in inequalities (77)-(79) remain strict, however, 
because they were obtained from the strict inequality (76). Thus all that is 
necessary is to set E m - ln = 0 everywhere following inequality (75). To see that 
the discussion from Eq. (72) onwards applies equally well to the stochastic pro- 
cess defined in Sec. 4, replace the numerator in Eq. (72) by trP t , where P f is 
the covariance matrix of that process, and then re-interpret the various other 
quantities following Eq. (72) accordingly. 


10 Conclusions 

The problem of finding a system of approximate evolution equations for the 
mean and covariance matrix of second-order stochastic processes defined by un- 
forced, nonlinear, conservative systems of ordinary differential equations with 
random initial conditions has been examined in this article from the viewpoint of 
energetics. A brief treatment of the stochastic initial-value problem for conser- 
vative nonlinear systems of ordinary differential equations was given, and it was 
used to show that the mean and covariance matrix of the resulting stochastic 
process are dynamically linked through an energy relationship. 

The second-moment closure equations are a nonlinearly coupled system of 
approximate evolution equations for the mean and covariance matrix of this 
process. An existence and uniqueness theory was given for these equations, 
based largely on existence and uniqueness theory for the stochastic initial- value 
problem. This theory was then used to show that, under appropriate hypotheses, 
the mean and covariance matrix whose evolution is given by the second-moment 
closure equations are the mean and covariance matrix of an additional, well- 
defined stochastic process. It was shown further that the second-moment closure 
equations are energetically consistent: the mean and covariance matrix whose 
evolution they define are dynamically linked through precisely the same energy 
relationship as that of the mean and covariance matrix of the original stochastic 
process. 

Several implications followed from this energetic consistency. One is that 
the total variance V) of each of the two covariance matrices, that of the orig- 
inal stochastic process and that of the one whose evolution is given by the 
second-moment closure equations, is strictly and identically bounded in time 
t. It was shown further that energetic consistency implies simple, identical, 
time- independent bounds on the ratio Vt/Vt 0 . Also it was shown that when 
the original conservative system is genuinely nonlinear, as defined in this arti- 
cle, total variance for the second-moment closure equations may be increasing, 
decreasing or stationary at times. 

Essential to the results of this article is that no assumption was made on the 
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initial probability distribution, apart from the existence of two moments. It was 
also argued that the normal distribution is not appropriate in general for conser- 
vative dynamics, because it requires the existence of realizations, with nonzero 
probability, having total energy larger than any given amount. Furthermore, 
for atmospheric and ocean dynamics there are mass-like and/or temperature- 
like state variables that are constrained by the dynamics to be positive. With 
this motivation, an hypothesis was introduced that requires the realizations to 
have bounded energy, with probability one, and appropriate state spaces were 
introduced to handle state variables that are bounded from below. Many of the 
results were illustrated with examples. 

The behavior of solutions of the second-moment closure equations was con- 
trasted with the behavior of solutions of the approximate system obtained by 
neglecting the nonlinear coupling term in the mean equation. This approximate 
system, usually derived under an assumption of normality, is at the heart of 
many current large-scale computational applications in atmospheric and ocean 
dynamics. It was shown that this approximate system is not energetically con- 
sistent, because the nonlinear coupling term is crucial for energetic consistency. 

The results of this article have left a number of open questions. For instance, 
it will be important to establish hypotheses under which solutions of the second- 
moment closure equations exist over arbitrarily long time intervals, in case solu- 
tions of the nonlinear dynamical system from which they are derived also have 
this property. It will also be important to develop efficient computational al- 
gorithms for implementing the second-moment closure equations in large-scale 
oceanic and atmospheric applications. Addressing these issues successfully will 
require the continued strong collaboration amongst physical scientists, compu- 
tational scientists and mathematicians that has been so fruitful in recent years, 
as this volume attests. 
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